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ABSTRACT 

Context. The interpretation that bimodal colour distributions of globular clusters (GCs) reflect bimodal metallicity 
distributions has been challenged. Non-linearities in the colour to metallicity conversions caused for example by the 
horizontal branch (HB) stars may be responsible for transforming a unimodal metallicity distribution into a bimodal 
(optical) colour distribution. 

Aims. We study optical/near-infrared (NIR) colour distributions of the GC systems in 14 E/S0 galaxies. 
Methods. We test whether the bimodal feature, generally present in optical colour distributions, remains in the opti- 
cal/NIR ones. The latter colour combination is a better metallicity proxy than the former. We use KMM and GMM 
tests to quantify the probability that different colour distributions are better described by a bimodal, as opposed to a 
unimodal distribution. 

Results. We find that double-peaked colour distributions are more commonly seen in optical than in optical/NIR colours. 
For some of the galaxies where the optical (g — z) distribution is clearly bimodal, a bimodal distribution is not preferred 
over a unimodal one at a statistically significant level for the (g — K) and (z — K) distributions. The two most cluster- 
rich galaxies in our sample, NGC 4486 and NGC 4649, show some interesting differences. The (g — K) distribution of 
NGC 4649 is better described by a bimodal distribution, while this is true for the (g — K) distribution of NGC 4486 GCs 
only if restricted to a brighter sub-sample with small K-band errors (< 0.05 mag). Formally, the K-band photometric 
errors cannot be responsible for blurring bimodal metallicity distributions to unimodal (g — K) colour distributions. 
However, simulations including the extra scatter in the colour-colour diagrams (not fully accounted for in the photomet- 
ric errors) show that such scatter may contribute to the disappearance of bimodality in (g — K) for the full NGC 4486 
sample. For the less cluster-rich galaxies results are inconclusive due to poorer statistics. 

Conclusions. A bimodal optical colour distribution is not necessarily an indication of an underlying bimodal metallicity 
distribution. Horizontal branch morphology may play an important role in shaping some of the optical GC colour 
distributions. However, we find tentative evidence that the (g — K) colour distributions remain bimodal in the two 
cluster-rich galaxies in our sample (NGC 4486 and NGC 4649) when restricted to clusters with small K-band photo- 
metric errors. This bimodality becomes less pronounced when including objects with larger errors, or for the (z — K) 
colour distributions. Deeper observations of large numbers of GCs will be required to reach more secure conclusions. 
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1. Introduction 

Globular cluster (GC) systems exhibit a bimodal optical 
colour distribution in the majority of luminous and interme- 
diate luminosity early-type galaxies (Ashman & Zcpf 1993, 
Elson & Santiago 1996, Peng et al. 2006). This is widely 
interpreted as being due to the presence of two old sub- 
populations 10 Gyrs) that differ in metallicity (Brodie 
& Strader (2006) and references therein). Metallicity bi- 
modality is obvious for the Milky Way (Zinn 1985, Bica et 
al. 2006) but seems less evident for our spiral neighbour 
M31 (e.g. Caldwell et al. 2011). 
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The presence of two peaks in the GC colour distribu- 
tions has been often taken as an argument in favour of the 
existence of two major epochs/mechanisms of star forma- 
tion in the host galaxies. There has been much debate in 
the literature over the past decades as to what is responsi- 
ble for this bimodality (see Brodie & Strader 2006). Among 
the scenarios for GC formation that account for it, one can 
list Ashman & Zepf (1992), Forbes et al. (1997), Cote et 
al. (1998), Beasley et al. (2002), Strader et al. (2005) and 
Rhode et al. (2005). All of the above assumed different for- 
mation channels for blue (metal-poor) and red (metal-rich) 
clusters. However, the properties of these cluster popula- 
tions are generally not too distinct. For instance, Peng et 
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al. (2008) studying the specific frequencies of 100 early-type 
galaxies from the Virgo Cluster Survey find basically the 
same trends when separating the GC systems in red and 
blue. Peng et al. (2006) analysing the GC systems colour 
distributions of the same data set find that there is great 
similarity between the GC peak-metallicity galaxy-mass re- 
lations for the two populations. This implies that the con- 
ditions of GC formation for metal-poor and metal-rich GCs 
could not have been too different. 

The complexity of the formation histories of early-type 
galaxies, within the hierarchical merging framework (e.g. 
Renzini 2006), are not expected to naturally produce the 
near universality of bimodal GC metallicity distributions. 
This is actually the case for models of GC formation built 
upon this hierarchical framework (e.g. Beasley et al. 2002) 
where bimodality occurs only after introducing a mech- 
anism that artificially truncates the formation of metal- 
poor GCs. Nonetheless, more recently, Muratov & Gncdin 
(2010) introduce a scenario in which the formation of the 
two metallicity sub-populations of GCs may be a natural 
outcome of the hierarchical theory of galaxy formation in 
some, although not the entire range of model realisations. 
Muratov & Gnedin (2010) prescribe the formation of GCs 
semi-analytically using assembly histories from cosmolog- 
ical simulations combined with observed scaling relations 
for the amount of cold gas available for star formation. 

The interpretation that colour distributions translate 
directly into metallicity distributions was challenged by 
Yoon et al. (2006). They demonstrate that non- linear 
colour-metallicity relations caused by the horizontal branch 
morphology (HB) may transform a unimodal metallicity 
distribution into a bimodal optical colour distribution. This 
issue has been investigated in more detail by Cantiello & 
Blakeslee (2007) , who conclude that combinations of opti- 
cal and near-infrared (NIR) colours are much less sensitive 
to this effect. Richtler (2006) had already shown that a 
flat metallicity distribution can result in a bimodal colour 
distribution using the Washington photometric system (ac- 
knowledging discussions with Boris Dirsch). 

A few observational attempts to address whether op- 
tical colour bimodality is really representative of metallic- 
ity bimodality exist in the literature, providing sometimes 
ambiguous or conflicting results. For example, Strader et 
al. (2007) find clear evidence for two metallicity sub- 
populations in the spectroscopic sample of 47 NGC4472 
GCs from Cohen et al. (2003). For the same set of data, 
though, the analysis of Cohen et al. do strongly favour 
bimodality. Kundu & Zcpf (2007) find an optical/NIR 
(/ — H) bimodal distribution for NGC 4486 in a small sam- 
ple of GCs from one NICMOS/HST pointing combined 
with HST/WFPC2 data. Moreover, Spitler et al. (2008) 
presented an optical/mid-infrared analysis using Spitzer 
Space Telescope for NGC 4594 and NGC 5128. Both galax- 
ies present a clear optical bimodality. While NGC 5128 
presents an obvious R-[3.6] bimodal distribution, more com- 
patible with the multimodal peaks that the spectroscopic 
work of Beasley et al. (2008) find (see also Woodley et 
al. 2010), NGC 4594 has a less clear bimodal distribu- 
tion in this baseline. In contrast, Alves-Brito et al. (2011) 
present spectra for over 200 GCs in the Sa NGC 4594 and 
find a clear bimodal distribution. More recently, Foster 
et al. (2010) and Foster et al. (2011) find very similar 
Calcium triplet values for red and blue GCs in NGC 1407 
and NGC 4494 despite their large colour difference. Since 



Calcium triplet is a metallicity indicator, similar values for 
red and blue GCs indicate similar metallicities. One possi- 
ble explanation given by the authors is the non-linear con- 
version between colour and metallicity. 

Blakeslee et al. (2010) simulate GC populations with 
both a mass- metallicity relation (blue-tilt) and a non-linear 
colour-metallicity relation and find bimodal colour distri- 
butions with a blue-tilt even though the metallicity distri- 
bution appears unimodal. 

Finally, the very recent works by Yoon et al. (2011a) and 
Yoon et al. (2011b) show that colours such as (u — g) and 
(u — z) have significantly less inflected colour metallicity 
relations than (g — z). They also show that the metallic- 
ity distributions obtained from inflected colour metallicity 
relations are strongly-peaked, unimodal and with a broad 
metal-poor tail, similar to that of the resolved field stars in 
nearby elliptical galaxies and those produced by chemical 
evolution models of galaxies (e.g. Bird et al. 2010). 

The purpose of this work is to investigate the nature 
of optical/NIR colour distributions of different GC systems 
in several early-type galaxies. We aim to shed light on the 
true nature of the metallicity distributions of GC systems. 

2. Observations and data 

The observations and data reduction techniques of the data 
used in this study are described thoroughly in Chies-Santos 
et al. (2011a). Here we briefly summarise the procedures 
applied. A sample of 14 early-type galaxies was imaged 
in the Ks-band (from now on referred to as only K) with 
LIRIS at the WHT and combined with archival ACS/HST 
F475W (~ g) and F814LP (~ z) images. The galaxies have 
Mb < — 19 and (m — Al) < 32. ACS images were reduced 
with MULTIDRIZZLE (Koekemoer et al. 2002) and LIRIS 
images with LIRISDR, in addition to standard IRAF rou- 
tines. 

The GCs were detected and had sizes (R e ff) measured 
in the ACS images with DAOFIND and ISHAPE (Larsen 
1999), respectively. Aperture photometry was performed in 
the g, z and K bands with PHOT (Stetson 1987). After the 
following criteria were applied to automatically detected 
sources, g < 23, 0.5 < (g - z) < 2.0, 1 < R cf t(pc) < 15 a 
careful visual inspection was performed where obvious non- 
cluster objects were removed. Also, sources that were too 
close together in the ACS images and that appeared as one 
bigger source in the LIRIS images, due to its lower reso- 
lution compared to ACS, were flagged. Finally, we caution 
the reader that the data of NGC 4382 and NGC 4473 were 
taken at non-photometric conditions. 

3. Integrated colours and horizontal branch 
morphology 

In this section we study the effect of the HB in the inte- 
grated colours of GCs. The colour of the HB varies abruptly 
between [Fe/H]= -0.6 and -0.9. At this [Fe/H] range the 
HB departs from the red-clump position (Lee et al. 1994). 
This is argued as the cause of the non-linear behaviour of 
the colour-metallicity relation (Yoon et al. 2006). Peng et 
al. (2006) with data from the Milky Way, NGC 4486 and 
NGC 4472 showed that the empirical transformation from 
metallicity to the optical (g — z) colour is clearly non- linear. 
Dirsch et al. (2003) had already called attention to the fact 
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that the colour-metallicity relation is non-linear and that 
starting from the bluest colours equidistant colour intervals 
are projected onto progressively larger metallicity intervals. 

In Fig. 1 (g — z), (jg — K) and (z — K) colour-metallicity 
relations are shown. For those relations we used 14 Gyr 
SPoT 1 (Raimondo et al. 2005, Biscardi et al. 2008) and 
YEPS 2 models. The models for (g — z) show a clear de- 
parture from linearity around [Fe/H]~ —0.5; the wiggly 
feature noticed by Yoon et al. (2006). This wiggle that is 
also visible in other commonly used colours such as (V — I) 
and (B — I) (Cantiello & Blakeslee 2007) is generated by 
the transition from blue to red HBs. This feature is pos- 
sibly also recognizable in (g — K), although very mildly. 
Such feature is not visible when considering (z — K). This 
last colour does however show a non-linearity at interme- 
diate metallicities more pronounced for the SPoT models 
with respect to the YEPS. This is probably due to an in- 
terpolation effect as the point responsible for this, [Fe/H] 
~ —0.4, is right at the transition between blue and red 
HBs. Nevertheless the sign of the second derivative of the 
(z — K) — [Fe/H] relation at the non-linearity is opposite 
from the one for (g — z) — [Fe/H]. This indicates the differ- 
ent nature of the non-linearities in these 2 colour spaces. 
If one would remove the [Fe/H] ~ —0.4 point, the (z — K) 
model would bear linearity The SPoT and YEPS models 
agree well for the (g — z) and (g — K) colour — metallicity 
relations. For the (z — K) colour — metallicity relation there 
is more disagreement between the models, with the YEPS 
models being systematically bluer in this colour than SPOT 
for a given metallicity. This also happens for (g — K) but 
to a smaller extent. 

In Fig. 2 the (g — K) — (z — K) diagram for the joint 
NGC 4486 and NGC 4649 GC sample is shown with differ- 
ent simple stellar population (SSP) models over plotted in 
different panels. Note the wavy feature the data presents 
around (g — K) ~ 3.2 and (z — K) ~ 2. A running median 
on the data points is over plotted in the top panel for bet- 
ter visualisation of the wiggly behaviour. The SPoT 14 Gyr 
SSP (middle panel) is over plotted where this non-linear 
feature is also present. SPoT models match the observed 
colour magnitude diagrams of star clusters fairly well over 
a wide range of ages and chemical compositions (Brocato et 
al. 2000, Raimondo et al. 2005). In fact, SPoT models are 
optimised to simulate the HB spread observed in galactic 
GCs and this may be one point of uncertainty. Sohn et al. 
(2006) has shown evidence that the most massive NGC 4486 
GCs may have different He abundances and thus different 
HB morphology than galactic GCs. The YEPS models are 
also shown in Fig. 2. As SPoT, the YEPS models consider 
the systematic variation in the mean colours of HB stars as 
a function of metallicity. Even though the SPoT and YEPS 
tracks in Fig. 2 do not fit the data well in the redder part 
of the diagram, they present a wavy feature approximately 
at the same location as the data: the location pinpointing 
the transition from blue to red HB morphology. Other SSP 
models, in which the detailed modelling of the metallicity 
dependent HB morphology of the SPoT/YEPS models is 
not included, do not show this behaviour (bottom panels 



SPoT models can be downloaded from www.oa- 
teramo. inaf . it /spot . 

2 YEPS models can be downloaded from 
http://web.yonsei.ac.kr/cosmic/data/YEPS.htm 
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Fig. 1. 14 Gyr SPoT and YEPS (g - z), (g - K) and (z - K) 
colour-metallicity relations according to the legend in the middle 
panel. The black short-dashed-line in the first panel corresponds 
to the empirical relation to the Peng et al. (2006) data for Milky 
Way and Virgo Cluster GCs (Blakeslee et al. 2010). 



for Padova08 3 and for Maraston 2005). The SPoT/YEPS 
models for the other 2-colour combinations also present this 
non-linear behavior, although it is not possible to see it in 
the data using (g — K) vs. (g — z) and (g — z) vs. (z — K) 
as clearly as it can be seen using (g — K) vs. (z — K). We 
suspect that this could be an effect of scatter as the K err 
is dominant. The errors in (g — K) and (z — K) are not ac- 
tually independent, but both dominated by the K err . This 
means that the uncertainties cause the data to scatter more 
or less diagonally in the plot. Moreover this diagram has a 
more similar scale than with the other combinations, where 
only one of the colours contains K. We refer the reader to 

3 Padova SSPs retrieved from the CMD 2.2 input form 
(http://stev.oapd.inaf.it/cmd), with Marigo et al. (2008) 
isochrones 
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Fig. 2. (z-K) vs. (g-K) for GCs ofNGC4486 and NGC4649, 
indicated by different symbols and with error bars in the bottom 
panel. Note the wavy feature the data presents around (z — K) ~ 
2 and (g — K) ~ 3.2. Upper panel: a running median on the data 
points; Middle panel: 14Gyr SpOT-Teramo and YEPS SSPs 
with a realistic treatment of the HB morphology (Raimondo et 
al. 2005, Yoon et al. 2006, Yoon et al. 2011a, Yoon et al. 2011b); 
Bottom panel: Maraston (2005) with blue and red HB at the 
lowest metallicities and Padova08 14 Gyr SSPs. See Fig. 16 of 
Chies-Santos et al. (2011b) for a (g — K) vs. (g — z) plot. 



Fig. 16 of Chies-Santos et al. (2011b) for a plot of (g - K) 
vs. (g- z). 

By excluding the [Fe/H] ~ -0.4 point for the SPoT 
models, the (z — K )-[Fe/H] would be more linear. However, 
a similar behaviour is found also for the YEPS models at 
approximately this metallicity. The (g — K) colour presents 
a near linear relation to metallicity and also has a broader 
baseline than (z — K). 



4. Colour distributions and bimodality tests 

In this section we investigate the behaviour of the (g — z), 
(g—K) and (z—K) colour distributions. They are shown for 
the different GC systems in Fig. 3. Note that while a great 
number of systems look bimodal in (g — z) they appear less 
and less bimodal as one moves on to the colours that should 
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Table 1. KMM outputs for the colour distributions of the galaxies: (1) 
galaxy and its colour distributions, (2) number of clusters assigned to the 
blue peak (-/V5), (3) number of clusters assigned to the red peak (N r ), (4) 
the mean of the blue peak (^b), (5) the mean of the red peak (/i r ), (6) 
the common width of the peaks (cr) and (7) the probability for rejecting 
a unimodal distribution (P(KMM)). 



in principle sample less the horizontal branch ((g — K) and 
(z — K)). The bimodality feature is generally the least ap- 
parent in the (z — K) histograms. To quantify whether the 
distributions are better described by bimodal as opposed to 
unimodal models, the KMM (Ashman, Bird & Zepf 1994) 
algorithm was applied to the data. The bimodal (blue and 
red peaks and their sum) distribution estimates returned 
by the code are shown in Fig. 3. The outputs of KMM 
for the different colours and galaxies are listed in Table 1. 
This test utilises the likelihood ratio test (LRT) to estimate 
the probability that the distribution of a number of data 
values is better modelled as a sum of two Gaussian dis- 
tributions than a single Gaussian. However, it obeys a \ 2 
statistics only when the two modes have the same variance 
(homoscedastic case). P(KMM) indicates the probability of 
rejecting a unimodal distribution in favour of a bimodal 
one, with a low probability indicating that a bimodal dis- 
tribution is preferred. 

The GMM (Muratov & Gnedin 2010) code was also ap- 
plied to the data. GMM is a robust generalisation of the 
KMM test. Besides modelling a bimodal distribution with 
different mode variances (heteroscedastic case), it uses three 
different statistics: LRT(x 2 ), distance from the mean peaks 
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(DD) and kurtosis (hurt). Moreover, it makes use of boot- 
strapping to derive uncertainties for the estimators of these 
different statistics. Finally, GMM outputs the probability 
that the three statistics give for rejecting a unimodal fit in 
favour of a bimodal one: P(x 2 ), P(DD) and P(kurt). From 
now on, we refer to P(x 2 ) as P(GMM). The DD statis- 
tic is a measure of how meaningful the split between the 
two peaks is. The kurt statistic is a measure of peaked- 
ness of the distribution. A positive kurtosis corresponds to 
a sharply peaked distribution whereas a negative kurtosis 
corresponds to a flattened distribution. A negative kurtosis 
is a necessary but not sufficient condition for bimodality. 
The DD and kurt are additional checks to support the re- 
sults of the LRT since the likelihood function is very sensi- 
tive to outliers far from the centre of the distribution and 
may reject a genuine unimodal distribution. 

Before running KMM and GMM objects out of the 
ranges 1 < (g — K) < 5 and 0.5 < (z — K) < 4 were 
removed from the distributions. These are objects far from 
the centre of the distributions and are very likely to affect 
the results of the bimodality tests by e.g. creating artificial 
peaks which will be interpreted as separate populations by 
these statistical tests. For NGC 4486, there is no object in 
this range, while for NGC 4649 three objects drop off the 
distribution with these restrictions. The outputs of GMM 
for the different colours and galaxies are listed in Table 2. 
For the galaxy NGC 4570, due to the low number of clus- 
ters, the test did not converge, and these estimates are not 
shown. 
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Fig. 3. (g — z), (g — K) and(z — K) colour distribution for the galaxies. Over plotted are the bimodal (blue and red peaks in 
corresponding colours) distributions returned by KMM and the sum of the blue and red peaks (in green). 
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In Fig. 4, histograms of the different probabilities re- 
turned by KMM and GMM are shown. The (g—z) probabil- 
ities of these tests are concentrated towards rather than 1, 
indicating the probable rejection of the unimodal distribu- 
tion. In fact, ~ 80% of the (g-z) distributions of P(KMM) 
and - 70% of P(GMM) are consistent with < 0.05. The 
P(DD) and P(kurt) are < 0.1 for - 40% and ~ 46% of the 
systems, respectively. Therefore a bimodal distribution is 
favoured for approximately half of the cases when the opti- 
cal colour is used. In contrast, the (g—K) probability values 
are spread between and 1. Only ~ 15% of the cases have 
P(KMM) and P(GMM) < 0.05. The P(kurt) is < 0.1 for 
only ~ 8% of the systems, while none have P(DD) < 0.1. 
Like (g — K), the (z — K) probability values are also spread 
between and 1. Compared to the former colour, a slightly 
larger number of systems have P(KMM) and P(GMM) < 
0.05 (~ 21% and - 38%, respectively). The P(DD) is < 
0.1 for only ~ 8% of the systems, while none have P(kurt) 
< 0.1. Note however, that generally for (z — K) and for a 
few cases of (g — K) these tests assign very few objects for 
one peak compared to the other, while in (g — z) the num- 
ber assigned to the different peaks is much more equal. In 
summary, the higher probability values assigned by KMM 
and GMM for (g — K) and (z — K) attest that bimodality 
does get less evident in (g — K) and (z — K) compared to 

(g-z). 

5. The cluster rich galaxies and simulations of the 
colour distributions 

The (z — K) and (g — K) colours, where bimodality ap- 
pears generally less pronounced, are also the colours with 
the higher photometric uncertainties (Chies-Santos et al. 
2011a). In the present section we investigate whether pho- 
tometric errors might be the reason bimodality is less ap- 
parent in such colours. 

First we take the three richest GC systems (NGC 4486, 
NGC 4649 and NGC 4552) to see how the distributions look 
when samples with restricted photometric errors are consid- 
ered. The colour distributions for these cluster rich galaxies 
are shown in Fig. 5 for different ranges of photometric un- 
certainties: K err < 0.5 (whole sample), K err < 0.1 and 
K err < 0.05. The outputs of KMM and GMM are shown 
in Tabs. 3 and 4 respectively. When applying these restric- 
tions, a greater fraction of blue objects drops off the sample, 
as compared to red objects. This is due to the colour de- 
pendant scatter, as discussed in Chies-Santos et al. (2011b). 
For the (g — z) distribution, all probability values are sig- 
nificantly low, indicating bimodality both for samples with 
K err < 0.1 and K err < 0.05. For the other colours the 
situation changes and distinguishing a unimodal distribu- 
tion from a bimodal one through KMM and GMM becomes 
very difficult. It can happen that a probability value goes 
from ~ 1 to ~ with the restricted sample, e.g. P(GMM) 
of (g - K) for NGC 4486. Also, it happens that most ob- 
jects are assigned to one of the peaks in the (g — K) and 
(z - K) KMM estimates of the K err < 0.1 and K err < 0.5 
samples (e.g. for NGC 4649). This exercise does not lead 
to a clear conclusion regarding the presence of bimodal- 
ity in the optical/NIR colour distributions, mainly due to 
the low number statistics when running the tests in the 
samples restricted in photometric uncertainties. As gener- 
ally assumed, if the (g — z) distribution is clearly bimodal 
as it is for NGC 4486 one would expect to see bimodal 
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Fig. 4. (g — z) [top panel), (g — K)(middle panel) and (z — K) 
{bottom panel) probability values from KMM and GMM. 



(g — K) and (z — K) distributions in the K err < 0.05 sam- 
ple. Instead, the optical/NIR distributions look much less 
bimodal. However, note that bimodality is somewhat signif- 
icant in (g - K) for NGC 4486 for the K err < 0.05 sample, 
with P(KMM)= 0.035 and P(GMM)< 0.001. 

5.1. Simulations including photometric scatter 

Here, we test whether bimodality could be blurred in 
(g — K) and (z — K) due to the high photometric uncertain- 
ties in K. As shown in Blakeslee et al. (2010) the most sub- 
tle details in the colour metallicity relation matter for the 
final transformed distribution. Therefore, to avoid unneces- 
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sary complications we do not use SSP models to transform 
the colours but assume linear transformations. A bimodal 
distribution (in this case, (g — z)) should maintain its form 
if linearly mapped to other colour spaces. We find the rela- 
tions between (g — z) and (g — K), and between (g — z) and 
(z — K) by fitting the best straight line by a procedure that 
gives more weight to the objects with less observational er- 
rors in both the abscissa and the ordinate using the GC 
systems of NGC4486 and NGC4649, as in Chies-Santos et 
al. (2011b). The relations are given by: 

(g -z) = 0.465 * (g - K) - 0.349 (1) 

(g -z) = 0.746 *{z-K)- 0.387 (2) 



many of the outcomes. Moreover, for NGC 4552, KMM does 
not converge for many cases, probably because of the low 
number of clusters. Therefore, the only galaxies for which 
we can robustly analyse the optical/NIR colour distribu- 
tions are NGC 4486 and NGC 4649. 

This analysis shows that if unimodality is preferred over 
bimodality in (g — K) and (z — K) it is not an effect of 
photometric errors in NGC 4486 and NGC 4649. KMM and 
GMM tests show that their (g — z) distributions are gen- 
uinely bimodal. The (g — K) distribution is unimodal for 
NGC 4486 while bimodal for NGC 4649. For the reddest of 
the colours, (z — K), bimodality is not significant for any 
of the galaxies. 



We linearly transform the (g — z) distribution to (g — K) 
and (z — K) using the inverted relations (1) and (2). In the 
top panels of Fig. 6 the results of this linear transforma- 
tion are shown for NGC 4486. Recalling the result of Chics- 
Santos et al. (2011a) that the photometric scatter in K, as 
measured by PHOT, is underestimated by a factor of 2, we 
add randomly a realistic K dependent scatter to the (g—K) 
and (z — K) transformed distributions. The contributions 
of the g and z errors to these colours are negligible if com- 
pared to the K errors. In order to add the error associated 
to each If -band magnitude to the colours we first determine 
an analytic relation between K vs. K err . The middle left 
panel of Fig. 6 shows the observed K and K err , the best fit 
relation as well as its one sigma deviation. The analytic re- 
lation is of the form K err — N x 10 a * K , where N and a are 
constants that are different for different galaxies. A random 
number drawn from a gaussian with dispersion correspond- 
ing to the measured one sigma deviation of the K vs. K err 
relation is calculated. As in Chies-Santos et al. (2011a) it 
was shown that the photometric errors are underestimated 
by a factor of 2, we multiply this random number by two 
and add it to the analytic K err . Following this procedure 
we associate a JT-band error to each GC according to its 
magnitude. 

One outcome of the randomly sampled realistic scatter 
is shown in the middle right panel of Fig. 6 as a scatter plot 
of large black symbols. For comparison, the parent distri- 
bution of errors is also shown as a scatter plot of small 
red dots. In the bottom panels the resulting (g — K) and 
(z — K) distributions after the scatter is added are shown. 
A realistic scatter is sampled randomly and added to the 
linearly transformed distribution 100 times to generate 100 
unique distributions. In order to have a sufficient number of 
clusters for a robust analysis, we do this only for the most 
GC-rich galaxies that show obvious bimodality in (g — z) 
and a nearly equal number of clusters in the blue and red 
peaks: NGC 4486, NGC 4649, NGC 4552 and NGC 4621. 

KMM and GMM ere run on the outputs of these sim- 
ulations. In Fig. 7, histograms of the probability values re- 
turned by these statistical tests are shown. For NGC 4486 
and NGC 4649 a linearly transformed bimodal distribution 
in (g — z) remains bimodal both in (g — K) and (z — K) 
even when realistic photometric scatter is added. This is 
attested by the statistically significant values of P(KMM), 
P(GMM) and P(kurt) (< 0.05) for the great majority of the 
cases. Nonetheless, the separation between the peaks is not 
as significant, P(DD) > 0.05. However, for NGC 4552 and 
NGC 4621 the high probability values show that photomet- 
ric scatter could blur bimodality in (g — K) and (z — K) in 
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Table 3. Results of KMM for NGC4486, NGC4649 and NGC4552 for the colour distributions with GCs K err < 0.1 and 
K err < 0.05: (1) galaxy, (2) number of clusters assigned to the blue peak (Nb), (3) number of clusters assigned to the red peak 
(N r ), (4) the mean of the blue peak (fib), (5) the mean of the red peak (/x r ), (6) the common width of the peaks (a) and (7) the 
probability for rejecting a unimodal distribution (P(KMM)). 




0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 

(z-K) (z-K) ( Z -K) 



Fig. 5. NGC 4486, NGC 4649 and NGC 4552 colour distributions in (g — z), (g — k) and (z — fc). The open histogram are the colour 
distributions for all clusters that make the final sample, according to the criteria outlined is section 2. The hashed histograms 
show the same distributions when only the clusters with K err < 0.1 are left in the sample and the filled histogram when only the 
clusters with K err < 0.05 are present. 
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Fig. 6. Illustration for NGC 4486 on how the simulations with 
realistic photometric errors were performed. Top panels: the (g — 
K) and (z — K) colour distributions transformed linearly from 
(g—z) using relations (1) and (2) respectively. Middle panels: the 
data and the modelled photometric scatter (K err ) as a function 
of K and one outcome of the randomly sampled photometric 
scatter (K errran ) as a function of K. Bottom panels: the (g — K) 
and (z — K) simulations with K errrarl added. 
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Fig. 7. Histograms of the probability values returned by KMM 
and GMM on the output of the simulations for the cluster rich 
galaxies with obvious (g—z) bimodality and nearly equal number 
of clusters assigned for the blue and red peaks; for (g — K) (left 
panels) and (z — K) (right panels). 



5.2. Simulations accounting for scatter in the colour-colour 
diagrams 

The analysis performed above is intended to reproduce 
the observed K-K err parameter space, and showed that 
the photometric uncertainty does not play a major role in 
"smoothing" the observed (g — K) and (z — K) colour distri- 
butions. However, it does not take into account the scatter 
in the (g — z) vs. (g — K) and (g — z) vs. (z — K) diagrams, 
that might be an intrinsic scatter, at least in part. 

In order to take into account the effect of the scatter in 
the colour-colour diagrams we performed two more sets of 
simulations for NGC 4486 and NGC 4649. For the second 
set of simulations, as before, we linearly transformed the 
(g — z) colour to (g — K) and (z — K) using Eqs (1) and (2). 
Following this, we calculated the residual for each GC from 
the transformation relations (1) and (2). For each galaxy 
we then simulated 100 colour distributions (both for (g — 
K) and (z — K)) by adding errors randomly pulled from 
the array of residuals with replacement to the transformed 
(g — K) and (z — K) disitributions. 

In the second paper of this series (see Fig. 14, Chies- 
Santos et al. 2011b) it is shown that the scatter estimated 
through a broken fit in the colour-colour diagrams (rather 
than a linear fit) is almost entirely consistent with the pho- 
tometric scatter. However, there we also note that there is 
still extra scatter found at bright magnitudes that may be 
attributed to age spreads or other factors. Following this, 
we performed a third set of simulations that consisted in 



using broken fit relations for transforming the (g — z) colour 
to (g — K) and (z — K). These are given by Eqs. (3), (4), 
(5) and (6). Eqs. (3) and (5) are valid for (g - z) < 1.187 
and Eqs. (4) and (6) for (g - z) > 1.187. 



(g - z) b = 0.260 *(g- K) b + 0.232 
( g - z) r = 0.340 * (g - K) r + 0.140 
(g - z) b = 0.292 *(z- K) b + 0.431 
(g - z) r = 0.437 *{z- K) r + 0.391 



(3) 
(4) 
(5) 
(0) 



Similarly to the first set, we ran KMM and GMM on 
the output of these set of simulations. Moreover, before 
running KMM and GMM we applied the same colour cuts 
(1 < (g - K ) < 5 and 0.5 < (z - K) < 4) as for the real 
data. In Fig. 8 we show the outputs of KMM and GMM 
of the third set of simulations, similarly to Fig. 7. For the 
third set of simulations these tests converged for > 96% 
of the cases for both KMM and GMM for both simulated 
colours. The (z — K) simulated distributions of NGC 4486 
converged for ~ 88% of the cases. Fig. 8 only shows the 
cases which converged. From the left panels of Fig. 8 it clear 
that there is a peak at 0.05 for P(KMM), P(GMM) and 
P(kurt). This indicated that for the (g — K) simulations 
of NGC 4486 and NGC 4649 bimodality is favoured over 
unimodaliy. Formally, ~44% of the simulations are found 
to be significantly bimodal considering P(KMM) for these 
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GC systems. If one considers P(GMM) these values are 38% 
for NGC4486 and 32% for NGC4649. 

However, the interpretation of the right panels of Fig. 8 
is slightly more tricky. First, the P(kurt) values are skewed 
to 1 especially for NGC4486, which attests the strongly 
peaked unimodal distributions. This, combined with tails of 
outliers and/or short baseline, will cause KMM and GMM 
to assign a high probability of bimodality to many of the 
simulated distributions. In fact, KMM (GMM) assigns a 
high probability of bimodality (< 0.05) for -52% (27%) 
and ~42%(17%) of the simulations respectively. KMM and 
GMM assign the great majority of the objects to one of the 
peaks and still give high probability of bimodality for sev- 
eral of the simulated (z — K) distributions. However, the 
parent (g — z) distributions show a nearly equal division 
between blue and red clusters. 

In summary, for the (g — K) distributions of NGC 4486 
and NGC 4649 it is reasonable to say that a bimodal (g — 
z) distribution when transformed to (g — K) will remain 
bimodal for a little less than half of the cases and at the 
same time reproduce the observed scatter in {g — z) vs. 
{g — K). However, one caveat to keep in mind is that the 
simulations based on the two colour diagrams are a worst- 
case scenario that assumes all the additional scatter comes 
from (g — K). The estimate as to how likely it is for a 
bimodal distribution to remain bimodal in the optical/NIR 
colours in the presence of a scatter that reproduces the 
observed spread in the colour-colour diagrams is uncertain 
for (z — K). We do not show the results of the second set 
of simulations as they are very similar to the third dataset. 

One might ask how well do the positions of the peaks of 
the (g — K) and (z — K) distributions agree with the (g — z) 
ones when "transformed back" to (g — z) using relations (1) 
to (6). For example, for NGC 4486 and NGC 4649, using the 
GMM values of Table 2 we find that the differences between 
the "transformed back" and the real distributions is very 
small, of the order of ~ 0.01 for relation (1), involving (g — 
K) and (g — z); but can be as high as ~ 0.2 for relation (2) 
involving (z — K) — (g — z). When the broken fit relations 
(3), (4), (5) and (6) are used the differences are < 0.06, 
being always higher for the relations involving (z — K). 

Given that (g — K) is a good metallicity proxy, if 
NGC 4486 has a unimodal (g — K) distribution, this anal- 
ysis suggests an underlying unimodal [Fe/H] distribution 
in the absence of no scatter in the colour-colour diagrams. 
This is contrary to NGC 4649, whose (g — K) colour and 
consequently its underlying [Fc/H], appear genuinely bi- 
modal. If colour bimodality is a consequence of metallicity 
bimodality as typically assumed, then one would have ex- 
pected more prominent optical/NIR bimodal distributions 
compared to the optical colour distributions, because the 
former trace metallicity better than the latter. However, bi- 
modality may be blurred by scatter and is detected in only 
half of the realizations. 

6. Discussion 

This study shows the appearance of the optical/NIR colour 
distributions of GC systems for over a dozen galaxies for the 
first time. The bimodality feature, clear in purely optical 
colours, is found to be less pronounced in redder colours. 
For the most GC-rich galaxies, NGC 4486 and NGC 4649, 
the K-band photometric scatter does not appear to be re- 
sponsible for the almost disappearance of this, once well- 
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Fig. 8. Histograms of the probability values returned by KMM 
and GMM on the output of the simulations for NGC 4486 (top 
panels) and NGC 4649 [bottom panels) for {g — K) and (z — K). 



established feature of GC systems. In the absence of any 
scatter in colour-colour diagrams such as (g — z) vs. (g — K) 
the non linearities caused by the HB morphology might in- 
deed be the cause of the optical colour bimodality as sug- 
gested by Yoon et al. (2006) at least in some of the GC 
systems. This seems to be the case for NGC 4486, but not 
for NGC 4649. However, before strong claims are made as 
to whether certain galaxies have or have not underlying 
bimodal distributions, data of better quality are necessary. 

In some galaxies, blue GCs tend to be redder at brighter 
magnitudes. This causes the observed blue-tilt (e.g. Strader 
et al. 2006, Mieske et al. 2006, Harris et al. 2006) fea- 
ture in the colour magnitude diagrams of certain GC sys- 
tems. One might argue that the blue-tilt combined with 
the higher photometric scatter present in the optical/NIR 
colours in comparison to the purely optical colour could be 
in part responsible for the disappearance of a clear bimodal- 
ity feature in optical/NIR colours of GC systems such as 
NGC 4486 (known to contain the blue-tilt). We argue that 
there is no reason to think so. By looking at the NGC 4486 
histograms in Fig. 5 for clusters with K err < 0.05, contain- 
ing only the brightest GCs of the sample, one can see that 
the (g — z) distribution shows a very clear bimodal feature 
while the optical/NIR distributions (especially (z — K)) do 
not show it as clearly. If the blue-tilt were to be (in part) 
responsible for the non-clear bimodal feature in the opti- 
cal/NIR colours it would have to also be in the optical 
one. Moreover, both blue-tilt and bimodality have been de- 
tected using (g — z) data. Furthermore this feature cannot 
be the main factor behind the absence of a clear (g — K) and 
(z - K) bimodality in NGC 4486 compared to NGC 4649. 
The GC systems of both NGC 4486 and NGC 4649 have a 
clear blue-tilt (Strader et al. 2006, Mieske et al. 2006). 

Cohen et al. (1998) presented a spectroscopic sample of 
GCs in NGC 4486. For a marginal preference of bimodality 
to be detected at the 89% significance level, they have to 
exclude a tail of very metal rich-GCs. Their [Fe/H] distri- 
bution is very narrow, in fact much more than that of the 
Milky Way, as shown in their Fig. 20. 

The fact that Kundu & Zepf (2007) found a clear (I-H) 
bimodal distribution for NGC 4486 contrasts with our re- 
sult for unimodality in (g — K). However, the (I — H) — 
metallicity relation, shown in Fig. 9 for SPoT and YEPS 
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models show a wavy feature that is responsible for project- 
ing equidistant metallicity intervals into larger colour bins. 
Could this be argued as the cause for the bimodal (I — H) 
distribution found by Kundu & Zepf (2007)? Note that the 
(I - H) dip in the colour distribution of NGC 4486 (Fig. 1 
of Kundu & Zepf 2007) occurs between 1.45— 1.65, exactly 
where the wavy feature predicted by the SPoT models is. 
As shown in Sect 3, the SPoT (g — if)-metallicity relation is 
far more linear in this regime. Moreover, the field of view of 
NICMOS, used in Kundu & Zepf (2007) is ~4 times smaller 
than that of LIRIS and ACS. Also, the number of clusters 
shown in the (I — H) distribution of Kundu & Zepf (2007) 
is roughly only 1 /4 of the number in this study and much 
more centrally concentrated, better sampling the red than 
the blue sub-population. Therefore the sample LIRIS /ACS 
sample is significantly different from the NICMOS one. 

To further investigate what could cause the bimodality 
seen in Kundu & Zepf (2007) we performed some simula- 
tions of unimodal metallicity distributions and transformed 
them to (I - H) using the 14Gyr SPoT model. This is 
done in order to show that the (/ — H) distribution might 
arise due to this specific colour-metallicity relation com- 
bined with the fact the sample is biased to metal-rich GCs. 
In the top panel of Fig. 10, examples of these simulations of 
unimodal metallicity distributions with different means for 
15000 (estimated NGC 4486 total number of GCs (Peng et 
al. 2006)) and 80 (number of GCs in the Kundu & Zepf 2007 
sample) GCs are shown. They are different gaussian distri- 
butions with the same dispersion = 0.5 and means = 0.2, 
0.7 and 1.2. The resulting respective colour distributions 
for these representative examples are shown in the bottom 
panel of Fig. 10. Note that the most metal-rich unimodal 
metallicity distribution becomes bimodal in (7 — H ) both 
for the 15000 and for the 80 GCs examples. Note also that 
the dips occur at the same (I — H) values (~ 1.5) of the 
distribution of Kundu & Zepf (2007). While the simula- 
tion that is more closely related to the NGC 4486 observed 
metallicity distribution (Cohen et al. 1998) is the interme- 
diate metallicity one, shown in black, the one shown in red 
matches better the observed data of Kundu & Zepf (2007) . 
This is probably because the (J — H) observed distribution 
is biased to more metal-rich values. It is well known that 
GCs in the centres of galaxies are more metal-rich than in 
their outskirts (e.g. Liu et al. 2011). One might argue that 
the (I—H) transformations shown in Fig. 10 are idealised as 
they do not contain photometric scatter. The uncertainty 
in (I - H) for the Kundu & Zepf (2007) study is found to 
be less than 0.1. Following this, a 0.1 gaussian scatter was 
added to the transformed (I — H) distributions. For the 
most metal-rich distributions (red histograms), bimodality 
continues to be seen in (I — H) for the great majority of 
the cases (with and without the addition of a 0.1 gaussian 
scatter). 

The near universality of GC colour bimodality, once well 
consolidated as direct implication for metallicity bimodality 
appears to go contrary to the evidence presented here. In 
the absence of any scatter in the colour-colour diagrams our 
analysis shows that the equivalence of bimodal colour and 
bimodal metallicity distribution is questionable and there is 
room for some speculation. Not all luminous galaxies neces- 
sarily have bimodal metallicity distributions. Each galaxy 
is particular and full of peculiarities, and its metallicity 
distribution should account for that. If in the Milky Way 
[Fe/H] bimodality is absolutely evident, this shows that in 




Fig. 9. SPoT and Yonsei 14Gyr (I - H) - metallicity relation. 
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Fig. 10. Simulations of unimodal metallicity distributions with 
different means for 15000 (upper left panel) and 80 GCs (bottom 
left panel) transformed to (I — H) according to the 14Gyr SPoT 
- (I — H) - metallicity relation. The resulting respective (/ — H) 
distributions are shown in the top and bottom right panels. 



our Galaxy the major star formation episodes occurred per- 
haps in less complicated ways. It might be that NGC 4649 
due to its obvious optical/NIR colour bimodality had a sim- 
ilar history to that of the Milky Way. CD galaxies, such as 
NGC 4486 and NGC 1399 (see the unimodal (I - H) distri- 
bution for NGC 1399 found by Blakeslee et al. 2012), are 
much more prone to interactions with their neighbours due 
to their physical position in the potential well of the galaxy 
cluster. Perhaps the quantity and/or strength of such in- 
teractions play an important role in determining the shape 
of the [Fe/H] distributions observed today. This picture of 
non- universality of [Fe/H] bimodality fits much better with 
the hierarchical merging paradigm. 

High S/N NIR/optical imaging in an 8-10m class tele- 
scope is of the utter most importance to further investigate 
this issue. Also, large spectroscopic datasets coupled with 
simulations for several GC systems would shed more light 
in the true nature of the metallicity distributions of GC 
systems in large galaxies. 
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7. Summary and conclusions 

We have analysed the (g — z), {g — K) and (z — K) colour 
distributions of GC systems in a sample of 14 E/SO galaxies. 
The results are summarised below: 

1. The data presents a non-linear feature around {z—K) ~ 
2 and {g — K) ~ 3.2, marking the transition from blue 
to red HB morphology. SPoT and YEPS models, with a 
realistic treatment of the HB also show a similar feature. 
According to these models for old ages the colour metal- 
licity relation does not present the prominent wiggle of 
Yoon et al. (2006) for the optical/NIR colours: (g - K) 
and {z — K). 

2. While the great majority of GC systems present an 
obvious bimodal distribution in (g — z), bimodality is 
clearly less pronounced in the optical/NIR distributions 
(g — K) and {z — K). The two most cluster rich galax- 
ies show some remarkable differences. While the GC 
system of NGC 4486 shows no obvious bimodality in 
{g — K) when all the GC sample is considered, the 
GC system of NGC 4649 does. However, if restricted 
to a brighter sub-sample with small K-band errors (< 
0.05 mag) the {g - K) distribution of NGC 4486 GCs is 
better described by a bimodal distribution. Simulations 
of the {g — K) and (z — K) distributions with realis- 
tic K-band errors suggest that the K-band errors can- 
not be the responsible for the blurring of genuinely bi- 
modal optical/NIR colour distributions in NGC 4486 
and NGC 4649. However, when taking into account the 
extra scatter present in colour-colour diagrams such as 
(g — K) vs. {g — z), we find that bimodality is indeed 
likely to be undetectable in over half of the cases for the 
(g — K) distributions of these two GC systems. 

3. The underlying metallicity distribution of the GC sys- 
tem of NGC 4649 appears to be a genuine case for bi- 
modality. However, for NGC 4486 the situation is less 
clear and bimodality is detected at a statistically sig- 
nificant level only for the brightest sub-sample of the 
clusters in {g — K) . This bimodality becomes less pro- 
nounced when including objects with larger errors, or 
for the {z — K) colour distributions. In the galaxy, cen- 
tre of the Virgo Cluster, the argument put forward by 
Yoon et al. (2006) might contribute for the clear optical 
colour bimodality. Higher S/N NIR imaging is strongly 
needed to understand whether the {g — K) lack of bi- 
modality for NGC 4486 could be due to scatter other 
than photometric. Also it would further constrain the 
nature of the optical/NIR colour distributions (and the 
underlying metallicity distributions) of the other, less 
GC-rich galaxies in the sample. 
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